#include "functions.h"
#include <random>
#include <chrono>

int main()
{
  typedef std::chrono::high_resolution_clock myclock;
  int n = 100;
  double** A = new double*[n];
  double** B = new double*[n];
  double** C = new double*[n];
  double** D = new double*[n];
  double** E = new double*[n];
  
  for(int i = 0; i < n; i++)
    {
      A[i] = new double[n]();
      B[i] = new double[n]();
      C[i] = new double[n]();
      D[i] = new double[n]();
      E[i] = new double[n]();
    }
  for(int i = 0; i < n; i++)
    {
      A[i][i] = B[i][i] = C[i][i] = D[i][i] = E[i][i] = 10.0;
      if(i == 0)
	{
	  A[i][i+1] = B[i][i+1] = C[i][i+1] = D[i][i+1] = E[i][i+1] = 1.0;
	}
	
      else if(i == n-1)
       	{
	  A[i][i-1] = B[i][i-1] = C[i][i-1] = D[i][i-1] = E[i][i-1] = 1.0;
	}
      else
	{
	  A[i][i+1] = B[i][i+1] = C[i][i+1] = D[i][i+1] = E[i][i+1] = 1.0;
	  B[i][i-1] = A[i][i-1] = C[i][i-1] = D[i][i-1] = E[i][i-1] = 1.0;
	}
    }
  //随机生成b，这里b暂取100以内的整数
  double* b = new double[n]();
  double* z = new double[n]();
  double* x = new double[n]();
  double* y = new double[n]();
  double* p = new double[n]();
  for(int i = 0; i < n; i++)
    {
      std::random_device e;
      b[i] = z[i] = x[i] = y[i] = p[i] = (double)(e()%100);
    }
  std::cout << "For the first matrix A, b is as following:" << std::endl;
  for(int i = 0; i < n; i++)
  std::cout << "b[" << i+1 << "] = " << b[i] << std::endl;

  //Root-Squaring Method
  myclock::time_point beginning_I1 = myclock::now();
  Root_Squaring(A, n);
  Lower_Triangular_Matrix(A, b, n);
  Upper_Triangular_Matrix(A, b, n);
  myclock::duration d_I1 = myclock::now() - beginning_I1;
  std::cout << "I.(1)Time of Root-Squaring Method:" << d_I1.count() << std::endl;
  std::cout << "I.(1)The Answer of Root-Squaring Method:" << std::endl;
  for(int i = 0; i < n; i++)
  std::cout << "x[" << i+1 << "] = " << b[i] << std::endl;
  for(int i = 0; i < n; i++)
    delete [] A[i];
  delete [] A;
  delete [] b;
  
  
  //Root-Squaring-pro Method
  myclock::time_point beginning_I2 = myclock::now();
  Root_Squaring_pro(B, n);
  Unit_Lower_Triangular_Matrix(B, z, n);
  Diag_Matrix(B, z, n);
  Unit_Upper_Triangular_Matrix(B, z, n);
  myclock::duration d_I2 = myclock::now() - beginning_I2;
  std::cout << "I.(2)Time of Root-Squaring-pro Method:" << d_I2.count() << std::endl;
  std::cout << "I.(2)The Answer of Root-Squaring-pro Method:" << std::endl;
  for(int i = 0; i < n; i++)
  std::cout << "x[" << i+1 << "] = " << z[i] << std::endl;
  for(int i = 0; i < n; i++)
    delete [] B[i];
  delete [] B;
  delete [] z;

  //Gauss-Elimination Method
  myclock::time_point beginning_I3 = myclock::now();
  Gaussian_Elimination(C, n);
  Unit_Lower_Triangular_Matrix(C, x, n);
  Upper_Triangular_Matrix(C, x, n);
  myclock::duration d_I3 = myclock::now() - beginning_I3;
  std::cout << "I.(3)Time of Gauss-Elimination Method:" << d_I3.count() << std::endl;
  std::cout << "I.(3)Answer of Gauss-Elimination:" << std::endl;
  for(int i = 0; i < n; i++)
    {
      std::cout << "x[" << i+1 << "]=" << x[i] << std::endl;
      }
  for(int i = 0; i < n; i++)
    {
      delete[] C[i];
    }
  delete [] C;
  delete [] x;

  //Column-Principle-Gauss-Elimination Method
  myclock::time_point beginning_I4 = myclock::now();
  Column_Principle_Gaussian_Elimination(D, y, n);
  Unit_Lower_Triangular_Matrix(D, y, n);
  Upper_Triangular_Matrix(D, y, n);
  myclock::duration d_I4 = myclock::now() - beginning_I4;
  std::cout << "I.(4)Time of Column-Principle-Gauss-Elimination Method:" << d_I4.count() << std::endl;
  std::cout << "I.(4)Answer of Column-Principle-Gauss-Elimination:" << std::endl;
   for(int i = 0; i < n; i++)
     {
       std::cout << "x[" << i+1 << "]=" << y[i] << std::endl;
       }
   for(int i = 0; i < n; i++)
     {
       delete[] D[i];
     }
   delete [] D;
   delete [] y;

   //LS-QR
   myclock::time_point beginning_I5 = myclock::now();
   LS_QR_Householder(E, p, n, n);
   myclock::duration d_I5 = myclock::now() - beginning_I5;
  std::cout << "I.(5)Time of LS-QR Method:" << d_I5.count() << std::endl;
  std::cout << "I.(5)Answer of LS-QR:" << std::endl;
  for(int i = 0; i < n; i++)
  std::cout << "x[" << i << "] = " << p[i] << std::endl;
  for(int i = 0; i < n; i++)
    {
      delete[] E[i];
    }
  delete [] E;
  delete [] p;
  
  n = 40;
  double** H = new double*[n];
  double** M = new double*[n];
  double** L = new double*[n];
  double** N = new double*[n];
  double** P = new double*[n];
  for(int i = 0; i < n; i++)
    {
      H[i] = new double[n]();
      M[i] = new double[n]();
      L[i] = new double[n]();
      N[i] = new double[n]();
      P[i] = new double[n]();
    }
  for(int i = 0; i < n; i++)
    {
      for(int j = 0; j < n; j++)
	{
	  H[i][j] = M[i][j] = L[i][j] = N[i][j] = P[i][j] = 1.0 / (i + j + 1.0);
	}
    }
  double* c = new double[n]();
  double* d = new double[n]();
  double* f = new double[n]();
  double* g = new double[n]();
  double* q = new double[n]();
  for(int i = 0; i < n; i++)
    {
      for(int j = 1; j < n; j++)
	c[i] = d[i] = f[i] = g[i] = q[i] = c[i] + 1.0 / (i + j);
    }

  //Root-Squaring Method
  myclock::time_point beginning_II1 = myclock::now();
  Root_Squaring(H, n);
  Lower_Triangular_Matrix(H, c, n);
  Upper_Triangular_Matrix(H, c, n);
  myclock::duration d_II1 = myclock::now() - beginning_II1;
  std::cout << "II.(1)[Hilbert-Matrix]Time of Root-Squaring Method:" << d_II1.count() << std::endl;
  std::cout << "II.(1)[Hilbert-Matrix]The Answer of Root-Squaring Method:" << std::endl;
  for(int i = 0; i < n; i++)
  std::cout << "x[" << i+1 << "] = " << c[i] << std::endl;
  for(int i = 0; i < n; i++)
    delete [] H[i];
  delete [] H;
  delete [] c;
  
  //Root-Squaring-pro Method
  myclock::time_point beginning_II2 = myclock::now();
  Root_Squaring_pro(M, n);
  Unit_Lower_Triangular_Matrix(M, d, n);
  Diag_Matrix(M, d, n);
  Unit_Upper_Triangular_Matrix(M, d, n);
  myclock::duration d_II2 = myclock::now() - beginning_II2;
   std::cout << "II.(2)[Hilbert-Matrix]Time of Root-Squaring-pro Method:" << d_II2.count() << std::endl;
  std::cout << "II.(2)[Hilbert-Matrix]The Answer of Root-Squaring-pro Method:" << std::endl;
  for(int i = 0; i < n; i++)
  std::cout << "x[" << i+1 << "] = " << d[i] << std::endl;
  for(int i = 0; i < n; i++)
    delete [] M[i];
  delete [] M;
  delete [] d;
  
  //Gauss-Elimination Method
  myclock::time_point beginning_II3 = myclock::now();
  Gaussian_Elimination(L, n);
  Unit_Lower_Triangular_Matrix(L, f, n);
  Upper_Triangular_Matrix(L, f, n);
  myclock::duration d_II3 = myclock::now() - beginning_II3;
  std::cout << "II.(3)[Hilbert-Matrix]Time of Gauss-Elimination Method:" << d_II3.count() << std::endl;
  std::cout << "II.(3)[Hilbert-Matrix]Answer of Gauss-Elimination:" << std::endl;
  for(int i = 0; i < n; i++)
    {
      std::cout << "x[" << i+1 << "]=" << f[i] << std::endl;
      }
  for(int i = 0; i < n; i++)
    {
      delete[] L[i];
    }
  delete [] L;
  delete [] f;

  //Column-Principle-Gauss-Elimination Method
  myclock::time_point beginning_II4 = myclock::now();
  Column_Principle_Gaussian_Elimination(N, g, n);
  Unit_Lower_Triangular_Matrix(N, g, n);
  Upper_Triangular_Matrix(N, g, n);
  myclock::duration d_II4 = myclock::now() - beginning_II4;
  std::cout << "II.(4)[Hilbert-Matrix]Time of Column-Principle-Gauss-Elimination Method:" << d_II4.count() << std::endl;
  std::cout << "II.(4)[Hilbert-Matrix]Answer of Column-Principle-Gauss-Elimination:" << std::endl;
   for(int i = 0; i < n; i++)
     {
       std::cout << "x[" << i+1 << "]=" << g[i] << std::endl;
       }
   for(int i = 0; i < n; i++)
     {
       delete[] N[i];
     }
   delete [] N;
   delete [] g;

   //LS_QR
   myclock::time_point beginning_II5 = myclock::now();
   LS_QR_Householder(P, q, n, n);
   myclock::duration d_II5 = myclock::now() - beginning_II5;
   std::cout << "II.(5)[Hilbert-Matrix]Time of LS-QR Method:" << d_II5.count() << std::endl;
   std::cout << "II.(5)[Hilbert-Matrix]Answer of LS-QR:" << std::endl;
   for(int i = 0; i < n; i++)
   std::cout << "x[" << i << "] = " << q[i] << std::endl;
  for(int i = 0; i < n; i++)
    {
      delete[] P[i];
    }
  delete [] P;
  delete [] q;
  
  return 0;
};
